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Abstract. We present a new numerical method for computing the pure-point 
spectrum associated with the hnear stability of coherent structures. In the 
context of the Evans function shooting and matching approach, all the relevant 
information is carried by the flow projected onto the underlying Grassmann 
manifold. We show how to numerically construct this projected flow in a 
stable and robust manner. In particular, the method avoids representation 
singularities by, in practice, choosing the best coordinate patch representation 
for the flow as it evolves. The method is analytic in the spectral parameter and 
of complexity bounded by the order of the spectral problem cubed. For large 
systems it represents a competitive method to those recently developed that are 
based on continuous orthogonalization. Wc demonstrate this by comparing the 
two methods in three applications: Boussincsq solitary waves, autocatalytic 
travelling waves and Ekman boundary layer. 



1. Introduction 

We introduce a new numerical method for solving high order linear spectral prob- 
lems by shooting and matching. The numerical construction of pure-point spectra 
is important in determining the linear stability of coherent structures. Examples 
of such structures are: ground and higher excited states of molecules in quantum 
chemistry (Johnson |56j, Hutson ^51j, Gray and Manopoulous |39j, Manopoulous 
and Gray [69] , Chou and Wyatt [23l [24] , Ledoux [64] , Ledoux, Van Daele and Van- 
den Berghe [65j . Ixaru [55j): nonlinear travelling fronts in reaction-diffusion such 
as autocatalysis or combustion (Billingham and Needham , Metcalf, Merkin and 
Scott !72], Doelman, Gardner and Kaper [3T], Terman [SB], Gubernov, Mercer, 
Sidhu and Weber [3T]); nerve impulses (Alexander, Gardner and Jones 2J), neural 
waves (Coombes and Owen ^25j): solitary waves or steady flows over compliant 
surfaces (Pego and Weinstein [82] , Alexander and Sachs [^ , Chang, Demekhin and 
Kopelevich [21], Kapitula and Sandstede [57], Bridges, Derks and Gottwald [14], 
Allen Allen and Bridges [5j); laser pulses (Swinton and Elgin [95^); nonlinear 
waves along elastic rods (Lafortune and Lega [B2j); ionization fronts (Derks, Ebert 
and Meulenbroek [28') or spiral waves (Sandstede and Scheel ^89]). 

For such problems the matching condition is typically a discriminant known as 
the Evans function (Evans [31], Alexander, Gardner and Jones [1]) or miss-distance 
function (Pryce [5^, Greenberg and Marietta [35]). It measures the degree of (possi- 
bly transversal) intersection of the stable and unstable solution subspaces satisfying 
the longitudinally separated far field boundary data. The stable subspace decays 
in the direction of wave propagation whilst the unstable subspace decays in the 
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opposite direction. Equivalently, the Evans function is the determinant of the set 
of solution vectors spanning both subspaces. With this end-goal matching condi- 
tion in mind, the problem boils down to how to numerically construct the solution 
subspaces in a robust fashion, as well as where to match longitudinally. This is espe- 
cially difficult for large scale problems. These might either emerge from high order 
systems, or more specifically, we envisage the stability of nonlinear travelling waves 
with multi-dimensional structure; for example wrinkled fronts travelling in a fixed 
longitudinal direction. The transverse structural information can be projected onto 
finite dimensional transverse basis, generating a large linear spectral problem posed 
on the one-dimensional longitudinal coefficient function set (see Ledoux, Malham, 
Niesen and Thiimmler [66j ) . Hitherto such large problems could not be solved by 
shooting and matching and were resolved by projecting the whole problem onto a 
finite basis and solving the resulting large algebraic eigenvalue problem. However 
recently, in the context of the Evans function, Humpherys and Zumbrun [49j pro- 
posed continuous orthogonalization as a viable approach to help make large scale 
problems amenable to shooting and matching. Here we provide our own answer. 

We propose the new Grassmann Gaussian elimination method (GGEM) which 
resolves several numerical problems all in one, in particular it: 

(1) Evolves the solution along the underlying Grassmann manifold avoiding 
representation singularities. 

(2) Retains analyticity in the spectral parameter. 

(3) Allows for matching anywhere in the longitudinal computational domain. 

(4) Has polynomial complexity, operations are of the order of the size of the 
system cubed. 

(5) Naturally evolves the solution in what in practice is the optimal coordinate 
representation (generated by optimal partial pivoting). 

For each property what is new, expected, proved, numerically observed, and its 
context? 

First, evolving the solution subspaces, considered as curves in the Grassmann 
manifold is not new for autonomous problems — see Hermann and Martin [44l |45l 
IMl mi SH], Martin and Hermann [TT] Brockett and Byrnes [18], Shayman [92], 
Rosenthal [57], Ravi, Rosenthal and Wang [55], Zehkin |100| . Abou-Kandil, Freil- 
ing, lonescu and Jank [1] and Bittanti, Laub and Willems [12]. Using Riccati sys- 
tems to solve nonautonomous spectral problems is also not new — see Johnson [56] , 
Hutson [51], Pryce [85], Manopoulous and Gray [69], Gray and Manopoulous [38] 
and Chou and Wyatt [23l [24] . Here Riccati systems correspond to the flow of the 
linear spectral problem projected onto the Grassmann manifold with a fixed coordi- 
nate patch representation — see Schneider [91] . Also see Schiff and Shnider [90] and 
Chou and Wyatt [24] who use this connection to integrate Riccati systems through 
singularities. 

That we have a numerical method that avoids representation singularities for 
nonautonomous systems is new. The idea is as follows. Given data that lies in 
a suitable coordinate patch of the Grassmann manifold, puUback to the Stiefel 
manifold. Evolve the solution one steplength along the Stiefel manifold either 
directly using a Runge-Kutta method or a Lie group method. Then project onto a 
suitable and possibly different coordinate patch of the Grassmann manifold using 
optimal Gaussian elimination. In this last step, the practical method we propose 
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picks a quasi-optimal coordinate patch in which to best represent the solution in 
the Grassmann manifold (see below and Section [5]). 

Note that on first inspection the Stiefel manifold is the direct natural setting for 
the stable and unstable solution subspaces. Afterall in each case we have to con- 
struct the full set of solutions to a large system of differential equations satisfying 
the correct respective asymptotic boundary conditions in the far field. Each solution 
set represents a curve in the Stiefel manifold of dimension commensurate with the 
size of the solution set. That the spectral problems are linear means that all the rel- 
evant spectral information can be reconstructed from the flow in the corresponding 
Grassmann manifold. That the matching condition is determinental, means that 
we only need the Grassmann flow information — see Martin and Hermann |71j and 
Brockett and Byrnes [18] where this reduction was first considered for autonomous 
linear control problems (in practice we will also need to retain a complex scalar 
field to ensure analytic dependence on parameters). This reduction is crucial be- 
cause long-range integration along the Stiefel manifold has been problematic (due 
to multiple distinct exponential growth and decay rates) and one of the simplest 
resolutions in the Evans function context was to use Pliicker coordinates — whilst 
ignoring the Pliicker relations (more on these below). 

Second, retaining analyticity away from the essential spectrum is standard for 
any shooting method; we prove analyticity in Section [S] This allows for a global 
search for eigenvalues in that region by numerically computing the change in ar- 
gument of the Evans function round any closed contour. Invoking the argument 
principle, this integer value represents the number of zeros of the Evans function, 
and hence the number of eigenvalues counting multiplicity, inside the contour; see 
Brockett and Byrnes [18J, Alexander, Gardner and Jones [2] and Ying and Katz 

Third, that our method allows matching anywhere in the longitudinal domain is 
new. We provide substantive numerical evidence. Previously, other than in trivial 
cases, most numerical practioners used the common-sense rule of thumb of integrat- 
ing the spectral problem from both ends of the longitudinal domain and matching 
at a point roughly centered on the front (which is assumed to be localized). When 
solving the linear problem with Pliicker coordinates, it was first important to rescale 
for the far field spatial behaviour to neutralize its total exponential growth. Inte- 
grating from the far field initial conditions (a subset of the spatial eigenvectors), 
the solution remains roughly constant until the coefficient matrix starts to reveal 
its nonautonmous character due to the integration step impinging on the front. Ac- 
curacy is retained whilst integrating through the front, but thereafter the problem 
becomes stiff. The issue is that the numerical methods cannot resolve the simul- 
taneous exponential growth and decay character associated with the other far-end 
stable and unstable subspaces. 

Fourth, having polynomial complexity is essential and we provide here a new 
alternate. After Humpherys and Zumbrun [IS] introduced their continuous orthog- 
onalization method in this context, which also has polynomial complexity, any new 
numerical spectral shooting method should have this property and also be competi- 
tive. Previous successful methods used Pliicker coordinates, ignoring the quadratic 
Pliicker relations (see Section [2]). They solved the flow for the corresponding linear 
vector field in the higher dimensional Pliicker embedding space. Details of this 
Pliicker coordinate or compound matrix approach can be found in, for example, 
Alexander and Sachs [3] , Brin [TBI HI] and Allen and Bridges [5] . Unfortunately 
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the number of Pliicker coordinates typically grows exponentially with the order of 
the original system, and so this approach cannot be used for medium to large or- 
der systems. However the continuous orthogonalization method of Humpherys and 
Zumbrun, and our method, are especially suited to large scale problems. 

Fifth, our method for dynamic practical optimal coordinate representation is 
new. Given data on the Stiefel manifold, for example generated by advancing the 
solution one steplength along the Stiefel manifold, how can we project down onto 
the Grassmann manifold using the best representation patch possible? The idea is 
as follows. The natural map projection from the Stiefel to the Grassmann manifolds 
is a linear fractional map (Milnor and StashcfF [74j Martin and Hermann JlJ). This 
map represents the action of equivalencing by transformations whose rank matches 
that of the Stiefel manifold (this takes us from the space of frames to the space of 
spaces spanning those frames). The Stiefel manifold has a non-square matrix rep- 
resentation. Projection onto the Grassmann manifold corresponds to equivalencing 
by a full rank submatrix of the non-square Stiefel matrix — this renders the corre- 
sponding submatrix as the identity matrix. We are free to choose which submatrix 
to equivalence by, each distinct choice corresponds to the matrix representation 
of a coordinate patch on the Grassmann manifold. We can use Gaussian elimina- 
tion, via elementary column operations, to render any given full rank submatrix of 
the Stiefel matrix as the identity matrix. The key is to try to pick the full rank 
submatrix which has the largest determinant — this corresponds to choosing the 
Grassmannian patch that gives the best representation for the projection from the 
Stiefel to Grassmann manifold. Ideally we would check the size of every full rank 
submatrix of the Stiefel matrix and equivalence by the one with the largest determi- 
nant. However this is an NP problem (equivalent to using the Pliicker coordinates 
described above). We provide a practical solution of polynomial complexity. The 
method maximises the pivot used at each step of the Gaussian elimination process. 
In the current context, we call it quasi-optimal Gaussian elimination. 

Our paper is organised as follows. In Section [21 we provide a tailored review of 
Grassmann manifolds and their representation. We then show, in Section [3] how 
flows generated by linear vector fields on the Stiefel manifold, produce a natu- 
ral flow on the underlying Grassmann manifold that is decoupled from the flow 
through the remaining fibres. We subsequently show how this leads to using Ric- 
cati systems to resolve spectra, but that singularities that develop in the Riccati 
flows present spectral matching problems. In Section[5]we introduce two new prac- 
tical approaches to avoiding these representation singularities. One is the idea 
behind our main method, the Grassmann Gaussian elimination method. The other 
is a modification of the Riccati approach that changes the coordinate patch when 
deemed necessary. Also in this section we show the connection between the Riccati 
and continuous orthogonalization approaches. We present our proposed Grassmann 
Gaussian elimination method fully in Section O including details of how in practice 
to choose the quasi-optimal Grassmannian coordinate representation patch. We re- 
view the Evans function in Section [6] and discuss further simple practical numerical 
refinements that retain analyticity and prevent potential numerical overflow. We 
then implement and compare all the competing numerical methods in Section [7| 
in three distinct applications. Finally in Section [5] we conclude and present future 
directions. 
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2. Review: Grassmann manifolds 

2.1. Grassmann and Stiefel manifolds. A fc-frame is a fc-tuple oi k < n linearly 
independent vectors in C". The Stiefel manifold Y{n,k) of fc- frames is the open 
subset of C"^'' of all fc- frames centred at the origin. The set of k dimensional 
subspaces of C" form a complex manifold Gr(n, k) called the Grassmann manifold 
of fc-planes in C" (see Steenrod [88j p. 35] or Griffiths and Harris [40l p. 193]). 

The fibre bundle tt: V(n, fc)— S'Gr(n, fc) is a principle fibre bundle. For each y in 
the base space Gr(n, fc), the inverse image Tr~^(y) is homeomorphic to the fibre 
space GL(fc) which is a Lie group — see Montgomery [TBI P- 151]. The projection 
map TT is the natural quotient map sending each fc-frame centered at the origin to 
the fc-plane it spans — see Milnor and Stasheff 74, p. 56]. 

2.2. Representation. Following the exposition in Griffiths and Harris [40], any 
fc-plane in C" can be represented by an n x fc matrix of rank fc, say Y G C"^*^. 
Any two such matrices Y and Y' represent the same fc-plane element of Gr(n, fc) if 
and only if Y' ~ Yu for some u S GL(fc) (the fc-dimensional subspace elements are 
invariant to rank fc closed transformations mapping fc-planes to fc-planes). 

Let i — {ii^...,ik\ C {!,..., n} denote a multi-index of cardinality fc. Let 
Yio C C" denote the {n — fc)-plane in C" spanned by the vectors {ej '■ j ^ i} and 

Vi = {Y e Gr(n, fc) : F n r,= = {0}} . 

In other words, Uj is the set of fc-planes Y e Gr(n, fc) such that the fc x fc submatrix 
of one, and hence any, matrix representation of Y is nonsingular (representing a 
coordinate patch labelled by i). 

Any element of U, has a unique matrix representation yio whose ith k x k sub- 
matrix is the identity matrix. For example, if i = {1, . . . , fc} then any element of 
can be uniquely represented by a matrix of the form 

/I ••• \ 

1 ••• 

1 

yk+l,k ' 
yk+2,k 

V yn,i yn,2 ■ ■ ■ yn,k / 

where j G C for i = fc + 1, . . . , n and j — 1^ . . . ,k. Conversely, a n x fc matrix of 
this form represents a fc-plane in Ui. Each coordinate patch Ui is an open, dense 
subset of Gr(n, fc) and the union of all such patches covers Gr(n, fc). For each i, 
there is a bijective map iy9i : Ui — > (£{n-k)k gjygj^ |-,y 

Vi'-Vi" ^y- 

Each ipi is thus a local coordinate chart for the coordinate patch Ui of Gr(n, fc). 
For all i, i', if y G Ui n Ui' and Mi_i' is the i'th fc x fc submatrix of ?/io, then 
2/(i')° ~ (wi,;')~^. Since Ui^i' represents the transformation between representative 
patchs and depends holomorphically on y\o , we deduce ipiOip~,^ is holomorphic. Note 
that Gr(n, fc) has a structure of a complex manifold (see Griffiths and Harris [40j 
p. 194]). Further the unitary group U(n) acts continuously and surjectively on 





ifk+l,! yk+1,2 
yk+2,1 yk+2,2 
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Gr(n, k). Hence Gr(n, k) is compact and connected. Lastly, the general linear group 
GL(n) acts transitively on Gr(n, k) and it is a homogeneous manifold isomorphic 
to GL(n)/GL(n - fc) x GL(fc) (see Chern [22, p. 65] or Warner J^, p. 130]). 

2.3. Pliicker embedding. There is a natural map, the Pliicker map, 

p: Gr(n, fc)^P(A'' C") 

that sends each /c-plane with basis Y = [Yi . . . Yfc] to Yi A . . . A Yfe; here P" de- 
notes the complex projective space of dimension n. See Griffiths and Harris |40j 
or Coskun [26] for more details. If we change the basis, the basis for the image 
changes by the determinant of the transformation matrix. Hence the map is a 
point in P(/\'^ C"). We can recover Y from its image Yi A ... A Yfc as the set of all 
vectors v such that wA Yi A . . . AY/c = 0. Further, a point of P(/\'' C") is in the image 
of p if and only if its representation as a linear combination of the basis elements 
of /\'' C", consisting of all possible distinct wedge products of a fc-dimensional ba- 
sis in C", is completely decomposable. Hence the image of p is a subvariety of 
P(/\'' C") of completely decomposable elements. It can also be realized as follows. 
A natural coordinatization of P(/\'^ C") is through the determinants of all the kx k 
submatrices of Y, normalized by a chosen minor characterized by an index i, hence 
P{/\'' C") = p(fe)^^. These minor determinants — the Pliicker coordinates — are not 
all independent, indeed, they satisfy quadratic relations known as the Pliicker re- 
lations (which may themselves not all be independent). The image of the Pliicker 
map p is thus the subspace of p(fe)^^ cut out by the quadratic Pliicker relations. 

3. Grassmannian flows 

3.1. Tangent space decomposition. Recall that we can consider the Stiefel man- 
ifold as a principle fibre bundle tt: V(n, fc) — > Gr(n, k). Our goal here is to charac- 
terize the induced decomposition of the tangent space TyV(n, fc) for Y e Y{n, fc). 
We can decompose the tangent space Ty V(n, fc) into horizontal and vertical sub- 
spaces (see, for example, Montgomery [THl P- 149]) 

TYV(n, fc) Hy e Vy . 

The horizontal subspace Hy is associated with the tangent space of the Grass- 
mannian base space, while the vertical subspace Vy is associated with the fibres 
homeomorphic to GL(fc). Let us choose the coordinate patch representation Ui 
for Gr(n, fc) for some i = {ii, . . . , ik} C {1, . . . , n}. Let Pi denote the projection 
matrix of size n x n that contains zeros everywhere except at positions {ii,ii) for 
/ = 1, . . . , fc where it contains ones. Note that one can additively decompose any 
given tangent vector V = PiV + PioV. Hence we have 

Hy = {PioV: y eTyV(n,fc)} C^'^^k)k ^ 

Vy = {PiV: y e TyV(n,fc)} = fll(fc) . 

3.2. Fibre bundle flow. Suppose we are given a vector field on the Stiefel manifold 

V{x,Y) = A{x,Y)Y, 

for any {x,Y) e M x Y{n,k) where A G Qi{n). Fixing a coordinate patch Ui for 
Gr(n, fc) for some i = {zi, . . . , ik} we can always decompose Y G V(n, fc) into 

Y = Viou, 
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where u £ GL(fc). Let a, b, c and d denote the i x i, i x i°, i° x i and i° x i° 
submatrices of A, respectively. 

Theorem 1. The flow governed by V{x,Y) generates a coupled flow in the base 
space Gr(n, fc) and fibres GL(fc); if Y = y^ou for a given fixed \, the flows in the 
coordinate chart variables y — (pi o y^o , and rank k transformations u, are 

y' = c + dy — y{a + by) and u' = {a + by)u, 

where we now think of a, b, c and d as functions of x, u and y. 

Proof. Using that Y = y\ou the ordinary differential system Y' = V{x, Y) becomes 

y[oU + y^ou' = (Ai + A^oy) u, 

where Ai represents the submatrix obtained by restricting the matrix A{x, yiou) to 
its ith columns. Applying the projections Pio and Pi to both sides of this equation, 
respectively generates the equations for y and u shown. Note that y is the projection 
of yio onto its i°th rows, as well as its image under the coordinate chart (fi. □ 

3.3. Riccati flow. A natural decoupling of the flow on the base space Gr(n, k) 
from the flow on the fibres GL(fc) occurs when the vector field V is linear, i.e. when 

V{x,Y) = A{x)Y, 

The following corollary is immediate from Theorem [1] 

Corollary 1. // the vector field V is linear so that A = A(x) only, then: 

(1) The flow on the base space Gr(n, k) decouples from the flow evolving through 
the fibres GL(fc) — the flow in the fibres is slaved to that in the base space; 

(2) For a fixed coordinate patch Ui index by i, in the coordinate chart variables 
y = ifi o y^o J the flow is governed by the Riccati differential system: 

y ^ c{x) + d{x)y - ya{x) - yb{x)y. 

Suppose we are required to determine the fiow generated by a linear nonau- 
tonomous vector field defined on the Stiefel manifold V(n, fc). The first conclusion 
in the corollary implies that all the relevant information is carried in the flow in the 
Grassmann manifold Gr(n, fc), and the flow through the fibres GL(fc) can be com- 
pletely determined a-posteriori from the Grassmannian fiow. The second conclusion 
suggests that if we fix a coordinate patch, then the flow in the Grassmannian can 
be determined from the solution to the Riccati system for y. If required, we can 
solve the differential system for u in Theorem[T]to determine Y , thereby completely 
resolving the flow generated by the linear vector field V on the Stiefel manifold. 
Provided y remains finite, this approach in fact works. 

The problem is that, though Y — yiou must be globally finite as it is generated 
by a linear vector field (with globally smooth coefficients), the Riccati solution y 
itself can become singular in a finite integration interval. Of course simultaneously 
the determinant of u itself becomes zero. The solution on the Grassmannian does 
not become singular. The issue is representation. (Note that the fiow on GL(fc) is 
linear but rank is not preserved because its coefficients depend on the Riccati flow.) 

The Riccati flow is a flow in a given fixed coordinate chart indexed by i, which 
is chosen at the start of integration. Given an initial element in V(n, fc), we pick a 
(good) coordinate patch U, for Gr(n, fc), this fixes the Grassmannian representation 
2/io . Projecting onto the i°th rows of j/io , or equivalently looking at the image under 
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the coordinate chart (y9i, generates y G ^(n-k)k^ rpj^^ Riccati flow is the flow in 
the Euchdean chart image space Each coordinate patch Ui is dense in 

Gr(n, fc). Therefore in numerical computations, the Riccati solution y is likely to 
leave and return to the patch Ui across any discrete integration step that staddles 
a representation singularity (generating a large but finite solution y the other side). 

With this in mind, Schiff and Shnider 90 suggested the following method that 
integrates through singularities in the Riccati flow. Fix a coordinate patch Ui with 
index i. The general linear group GL(?i) acts transitively on V(rt, k) (and also 
Gr(n, A:)): the left Lie group action A: GL(7i) x Y{n,k) Y{n,k) is defined by 
A: (5, Y)^ SY. For Yq G V(n, fc) we set 

Ay„:S^SYo. 

The Mobius Lie group action /i^,-, : GL(n)^C'^"^'^^'^ is defined by fiy^ : S (p\on\o 
^ip^^oyo ° where TTi : V(n, fc)^Ui is the quotient map tt; : Y y^o. Explicitly, if 
Siy represents the i x i' submatrix of S, we have: 

^J,yo■■ S (S'i.i + S'i,i°yo)(5'io4 + S'io,ioyo)~^ 

Thus, given data yo in the Euclidean chart image space 

^{n~k)k^ puUback to the Lie 
group GL(7i), to the identity element /„, using the Mobius Lie group action map 
fiyg. Advance the solution across one integration step in the Lie group generating 
the element S G GL(n). Schiff and Shnider used a Neumann/Runge-Kutta method 
to do this, but a Lie group method could also be used. Now push forward to 
y — o S* e £(ri-k)k ^gj^g Mobius Lie group action map. This takes you back 
to the chart corresponding to the original coordinate patch Ui . 

This method integrates through singularities in the Riccati flow. However, we 
are still left with another associated practical problem. If the linear vector field V 
depends on a parameter we wish to vary, the singularities in the Riccati flow can 
drift in the domain of integration and impinge on the matching position. 

4. Practical Grassmann integration 

Our goal in this section is construct numerical methods that integrate the flow 
associated with the push forward of the linear vector field V onto the Grassmannian 
manifold, whilst avoiding the representation singularities that occur in the Riccati 
approach. The solution is to change patch continuously in some optimal fashion or 
whenever the coordinate patch becomes a poor representation. 

The following diagram helps map out the two new strategies we suggest. 

^(„-fe)fe tl > y v(n, k) ^^"""^ ) GL(7i) — > Q\{n) 



GGEM 



Magnus 



^in-k)k ! U. ; Q^^^ V(n, k) GL(n) ; Qi{n) 



4.1. Continuous optimal patch evolution. The idea behind the Grassmann 
Gaussian elimination method (GGEM) is this. Given data in Gr(n, k) that lies 
in a given patch Yq £ Ui identified by i, puUback to the Stiefel manifold using 
the identity map — note Ui C V(n, fc) — so lo G Y{n,k). Advance the solution one 
integration step along the Stiefel manifold V(n, fc), using say a classical Runge- 
Kutta method, thus generating the element Y e V(n, k). 
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Now with the next step solution on the Stiefel manifold, we use quasi-optimal 
Gaussian elimination with partial pivoting (QOGE) to decompose Y — ?/(i/)ou and 
project onto the coordinate patch Uj' producing the element y(i')o e Ui'. The quasi- 
optimal Gaussian elimination process (described in Section picks out a suitable 
coordinate patch Uj' to represent the solution, which may be different than the 
original patch Ui. 

To ensure we numerically remain with the Stiefel manifold, rather than use a 
Runge-Kutta method, we might use a Lie group method as follows (see Munthe- 
Kaas [77] )■ Fullback the data Yq G V(n, k) from the Stiefel manifold to the general 
linear group GL(n), via the action map Ayp. The corresponding element in GL(n) 
is naturally the identity element /„. Subsequently pullback, via the exponential 
map, to the zero element in the corresponding Lie algebra, i.e. o € 0l{n). Evolve 
the solution on the Lie algebra to cr G 0i{n) using the Magnus expansion (Mag- 
nus [68]; also see Iserles, Munthe-Kaas, N0rsett and Zanna [54]). Pushforward 
from g[(n) to GL(n), via the exponential map, producing the Lie group element 
S = expcr S GL(n). Now pushfoward to the Stiefel manifold Y{n,k) via the Lie 
group action map Ay^ . For more details on Lie group methods for Stiefel manifolds, 
see Krogstad [61] and Celledoni and Owren [20] . 

4.2. Riccati flow with patch swapping. Since we are interested in constructing 
the flow generated by the push forward of the vector field V onto the Grassmann 
manifold, we can avoid singularities in any given Riccati system chart flow associ- 
ated with a given Grassmannian coordinate patch, by simply changing patch when 
the solution representation appears to become poor. In particular, we could either 
change the coordinate patch when the: 

• Norm \\y\\oD becomes too large (the easier and preferred approach we take); 

• Determinant det u becomes too small (this involves constructing u as inte- 
gration proceeds and we want to avoid carrying information unnecessarily) . 

Hence if say \\y\\oa exceeds a prescribed tolerance at the end of one step, to change 
patch, we apply the quasi-optimal Gaussian elimination method (QOGE) to the 
matrix yio = ip~^ o y. This identifies a new patch and index i' to use for the next 
set of successive steps until |l2/||oo becomes too large again, and so forth. 

4.3. Drury— Oja flow. There is a close connection between the Riccati flow in 
Corollary [1] and the continuous orthogonalization method of Humpherys and Zum- 
brun 49J, proved by direct comparison. 

Lemma 1. Ify £ ([^{■n-k)k g^UgjiQg ifi^ Riccati flow y' ~ c{x)+d{x)y~ya{x)—yb{x)y 
and u € U(fc) satisfies u' — —y^{c + dy)u, then for any index i, we have that 
Q = yi°u satisfies the Drury-Oja flow: Q' ~ {In — QQ^)AQ. 

Humpherys and Zumbrun [49j derive this flow by a Q-R-decomposition of the so- 
lution Y = QR G V(n, k) to the linear, spectral, globally bounded flow Y' = A{x) Y. 
We can think of this continuous orthogonalization method as generating an approx- 
imate flow on the Grassmann manifold whilst evolving the coordinatization (see 
Edelman, Arias and Smith [33j; Bindel, Demmel and Friedman [H]). It is also 
known as a Drury-Oja flow (see Drury [32], Oja [80], Yan, Helmke and Moore [98] . 
Bridges and Reich [15], Died and Van Vleck [30] and Hairer, Lubich and Wan- 
ner [421 P- 136]). The determinant of the upper triangular matrix, deti?, grows 
exponentially. Thus, since it is nonzero in the far field, we know it remains nonzero 
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in the whole integration interval K. Consequently, Q ~ YR^^ is globally finite on 
M, i.e. there are no singularities in the Drury-Oja flow. 

There is a natural su(n) Lie algebra action on the Stiefel manifold of orthonormal 
A:-frames, Vo(ti, fc), generated by the map (cr, Q) i— > exp((T) Q, for Q e Vo(rt, fc) and 
a e su(n). If u o Q = (/ — QQ^) A, the flow on the Lie algebra 5u(n) that generates 
the Drury-Oja flow on Yo{n,k) is governed by cr' = dexp~ o u(exp((T)(5o) ■ We 
could use this to construct a numerical method that preserves orthonormality of Q. 



5. Grassmann Gaussian elimination method 

5.1. Algorithm. The Grassmann Gaussian elimination method using a Lie group 
method on the Stiefel manifold (GGEM-LG), proceeds as follows: 

(1) Suppose initially we are given data j/jo (a;,„) in the coordinate patch Ui,„. 

(2) Across the integration interval [xm, Xm+i], compute cr^ using the Magnus 
expansion — we recommend the fourth order Magnus expansion. 

(3) Compute Y^+i = expert • j/jo (xm). 

(4) Apply quasi-optimal Gaussian elimination (QOGE) (outlined below) to 
Fm+i; this generates the solution yjo^^^ (xm+i) hi the coordinate chart 
^Xm+n ^^'^ the rank k matrix Um+i we have effectively equivalenced by. 

Across the integration interval [xm, Xm+i], we could generate Ym+i from yi^{xm) 
by solving the flow on the Stiefel manifold V(n, k) using a classical Runge-Kutta 
step (GGEM-RK). Both algorithms are summarized in the following diagram. 



GGEM 



id 



QOGE 



(Ayjo (.„))* 



RK 



Magnus 



cxp 



How much of Um+i do we retain at each step? We address this in Section [631 



5.2. Quasi-optimal Gaussian elimination. To project the endpoint solution 
Y„i+i in the Stiefel manifold onto a quasi-optimal coordinate patch, say Ut^^j , 
we use a Gauss- Jordan approach. This entails using elementary column operations 
and optimal pivoting, until a specified fc x /c submatrix of Y^^+i becomes the identity 
matrix (we naturally assume k > 2). To be more explicit about what we mean by 
optimal pivoting, we outline the procedure: 

(1) Look for the largest term in magnitude in Y'-^^ := Y^+i- Nominate this 
term as the pivot term and suppose this occurs in row zi, column ji. Use 
elementary column operations with this term, Pi-^j-^ := Y^pj_^, as the pivot 
to render all the other {k — 1) terms in that row equal to zero. Finally, 
use the elementary column operation of scalar multiplication to render the 
pivot term itself equal to one. Let us call the resulting n x k matrix Y^"^^; 
we can write (see for example Meyer [73]) 



y(2) ^ y(i) jjd)^ 
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where the elementary column operations performed are encoded in the el- 
ementary matrix U^^\ To be precise, we can write 

r/(l) _ ^(1) . . . ^(1) ^(1) ...^(1) ^(1) 

where for £ G {1, . . . , k}/{ji} the elementary matrix E^^^ encodes the fol- 
lowing elementary column operation on column ce: 

ce^ Ce . 

Piiji 

For i = ji the elementary matrix Ej^^ encodes the elementary column 
operation Cj^ Cj^/pi^j^. In summary, the matrix F^^^ has value one 
at position («i, ji) and otherwise zeros in row ii; and we have detC/^^^ = 
(ftiii)"^- 

(2) In the (n — 1) x (fc — 1) submatrix of F^^^ identified by excluding row ii 
and column ji, look for the largest term in magnitude. Again nominate 
this term as the pivot term and suppose this occurs in row 12 (which will 
be distinct from ii), and column j2 (which is distinct from ji). Here 12 
and j2 refer to the row and column relative to the original n x k matrix 
y'^-*. Use elementary column operations with this term, Pi^j^ := ^ij^]^' ^ 
the pivot to render the terms in row ^2 and columns {1, . . . ,k} / j^} of 
y*^^-* equal to zero. Again use the elementary column operation of scalar 
multiplication to render the pivot term itself equal to one. Let us call the 
resulting n x matrix Y^^^; we can write 

y(3) ^ y(2) jj(2) ^ y(l) jj(2)^ 

where the elementary column operations performed on Y^'^^ are encoded in 
the elementary matrix J7(2). Again, to be precise, we can write 

Tj{2) _ ^(2) . . . ^(2) ^(2) _ _ _ ^(2) ^(2) _ _ _ ^(2) ^(2) 

where for ^ e {1, . . . , ^'2} the elementary matrix Ef'^ encodes the 

following elementary column operation on column q: 

Cl^ Cl Cj^. 

Note that in our expression for t/*^^) above we could have either ji < 22 or 
ji > 22- For (. = j2 then encodes the elementary column operation 

— * Cj^/pi^j^. In summary, the resulting matrix y*-"*-* has value one at 
positions (*i, ji) and (12,^2) and otherwise zeros in row ii, and zeros in row 
12 in columns {1, . . . , 72}; and we have AetU'^'^^ = {Pi2j2)^^- 

(3) Continue this process. Focus on the (n — 2) x (fc — 2) submatrix of F^^^ 
identified by excluding the rows i\,i2 and columns j\ , j2 ; look for the largest 
term in magnitude. Nominate this term as the pivot term — suppose it 
occurs in row and column J3, relative to the original n x k matrix Y^^\ 
and so forth. On completing the final fcth step in this process, we will have 

y(fe) C/(l){7(2) ... [/W^ 
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where 

dct(;7(i) c/(2) . . . t/W) = . . •K,,,)-^ 

The final n x k matrix Y^''^ will have ones in positions {ie,je), for i = 
1, . . . ,k. In row ig it will have zeros in columns {1, . . . , A:} /{ji, ■ ■ ■ , je}. We 
set 

im+i := {«i,«2, • • • ,«fe}- 

(4) Perform column swaps in so that column ji becomes column 1. column 
j2 becomes column 2, and so forth so that finally column jk is forced to be 
column k. These column swaps can be encoded in the elementary matrix 
S with det S = (-l)#{s^aps} r^j^^ resulting matrix Y^''^ is given by 

y(fc) _ y(fe) Y,_ 

The im+i X {1, . . . , fc} submatrix of Y^'^^ given by 

L := [y* ^]i„+ix{i,...,fc} 
is lower triangular with ones on the diagonal. Finally we set 
2/i= := y('=) L-\ 

In practice of course, we do not compute L^^, but continue performing 
elementary column operations on y*^*^^ so that the submatrix L becomes the 
identity matrix, thus generating Vi^^^ ■ Hence we have effectively performed 
the decomposition 

Ym+l = • Um+1, 

where Um+i = L • • • (C/^^))"^ e GL(A;), and in particular 

uet Um+1 — \ -LJ Aiji PikJk ■ 

Note that if Ym+i G V(n, fc) then the entries in yi^^^ and Um+i: produced as 
a result of this process, will all be finite. Further if Ym+i depends analytically 
on a parameter, then the product yi<^ ^ ■ Um+i, naturally does as well. Lastly we 

remark that we could have performed alternative elementary column operations of 
the form ci — s- pijCg — Y^^^Cj (we do not include the scalar multipication operations 
here) with the result that det Um+i = (-1)#{~> {p^,nY-''{p^^nf-'' ■ ■ ■ ipikikf- 



5.3. Complexity. The complexity of the quasi-optimal Gaussian elimination al- 
gorithm, dominated by the search for the largest elements in the successively de- 
creasing submatrices of Ym+i, is of order nk"^. The method is a practical approach 
to maximize the determinant of the k x k submatrix removed from ym+i- It will 
not in general choose the submatrix with the largest determinant — hence the label 
quasi- optimal. This could be achieved by searching through all the k x k subma- 
trices of Ym+i, i.e. all the PKicker coordinates, but this has complexity of order n 
choose k. An interesting question here is whether there is an efficient way to use 
the Pliicker relations to reduce this complexity? 
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6. Spectral problems 

6.1. Linear Stiefel flow. Consider the linear spectral problem on M: 

Y' ^A{x-\)Y 

We assume there exists a subdomain C C containing the right-half complex 
plane, such that for A G O there exists exponential dichotomies on K~ and K"*" 
with the same Morse index k in each case (see Henry [l^ and Sandstede [88]). Let 
(x; A) € V(n, k) denote the matrix whose columns are solutions to the spectral 
problem and which span the unstable manifold section at x e [— oo,-t-oo). Let 
Y^{x]X) S Y{n,n — k) denote the matrix whose columns are the solutions which 
span the stable manifold section at a: e (— oo, +oo]. 

6.2. Matching. The values of spectral parameter A G $7 for which the columns 
of Y~ and columns of y+ are linearly dependent on R are pure-point eigenvalues. 
The Evans function D{X) is the measure of the degree linear dependence between 
the two basis sets Y~ and y+, i.e. of the degree of transversal intersection between 
the unstable and stable manifolds (see Alexander, Gardner and Jones 2j; Nii 79]): 

D{X)=e-^o^'^i^'^^<iidet{Y-{x;X) Y+{x;X)). 

It is analytic in f2. In practice we drop the non-zero, scalar exponential prefactor 
and evaluate the Evans function at a matching point x*. 

There are other matching criteria measuring the degree of intersection between 
subspaces that do not use the determinant: for example computing the angle be- 
tween subspaces as suggested by Bjorck and Golub [13] or computing the smallest 
eigenvalue as suggested by Hutson [51] and Ixaru [55] . Both these latter techniques 
might be important for large systems when computing the determinant could be an 
unstable process, indeed, we investigate them in this context in Ledoux, Malham, 
Niesen and Thiimmler [66]. However in both cases the magnitude of a function of 
the spectral parameter is computed, whose zeros correspond to eigenvalues. Hence 
we must search for touchdowns to zero in the complex parameter spectral plane 
which can be problematic. For the examples we consider here, which are not too 
large, using the determinant suffices. 

6.3. Initialization. We construct the nxk matrix Fq" (A) whose columns are the k 
eigenvectors of A{—oo] A) corresponding to eigenvalues with a positive real part (see 
Humpherys and Zumbrun [IS] and also Humpherys, Sandstede and Zumbrun |SD] 
for how to preserve analyticity with respect to the spectral parameter A). In prac- 
tice, integration starts at x = £- for some suitable, usually negative, value of i-. 
Analogously we construct the n x [n — k) matrix Y^{X) whose columns are the 
n — k eigenvectors of A{+oo; A) corresponding to eigenvalues with a negative real 
part. Again, in practice, we integrate backwards from x = £^ for some suitable, 
usually large and positive, value of £+. 

6.4. GGEM matching and analyticity. To compute Y^{x^,■,X) we start with 
Yi^(X) at X = i± and integrate centrally towards x — x^. Our goal in this section 
is to show that we only need the determinant of the rank k transformations in the 
GGEM method described at the beginning of Section [5l to retain analyticity for 
the Evans function in O. Since the procedure is the same in both intervals [^±,a;*] 
we will describe it for the generic interval [£, a;*] for Y{x; A) G V(n, k) starting with 
value Yo{£;X) a.t x = £. Suppose we use M successive computation subintervals 
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[xm,a;m+i] in [t,x*\ where Xm ■=(--{- m (x* — (.)/M. We label the nodal solution 
values at Ym{\). 

At the start xq = I, perforin quasi-optimal Gaussian elimination (QOGE) on 
lo(A) := 1^0 A) to obtain the decomposition 

yo(A) = t/io(a;o;A)[/o(A). 

As we shall see, we do not need to actually store t/o(A); but only det f/o(A). 

Let 'S'm,rri+i (A) denote an approximation to the flow-map across [xmi-^^m+i] to 
the linear system Y' = A{x\X)Y . We assume that Sjn,m+\{}^ preserves analytic 
dependency on A, so that the next step solution value Ym+ilA) = >S'm.m+i(A) ym(A) 
analytically depends on A if Ym{^ does. For example, in the case of GGEM-LG 
then Sjn,m+\{}^ = expcTTO and most straightforward Magnus based integrators will 
naturally preserve analyticity with respect to A. Similarly most simple Runge- 
Kutta methods used to generate S'm.m-i-iCA), or directly generate the next step 
solution value y^+i(A), will preserve analyticity. 

Our numerical procedure would proceed as follows. Across [xcxi] we have 

5o,i(A)yo(A) = (5o,i(A)yio(aro;A)) C/o(A) 
= t/ij(ari;A)C/i(A)C/o(A), 

where in the last step we applied QOGE to the n x k matrix S'o,i(A)2/io(a;o; A). 
Subsequently, across [a;i,a;2] we have 

5i,2(A) (yio(a;i;A)C/i(A)C/o(A)) = (5i,2(A) (xq; A)) Ui{X)Uo{X) 

= yio{x2;X)U2{X)Ui{X)Uo{X), 

where we applied QOGE to Si^2{X) yi°{xo; X). Repeating this argument across 
the subsequent intervals [xrmXm+i] for m = 2, ...,M — 1, we get the following 
approximation y(a;»; A) to Y{x^,; A): 

y(a;,;A) =^m-i,m(A) • • • ^oa(A) yo(A) 
= yi^(a;M; A) ;7m(A) ••• Uo{X) 
= yio^{xM;X)U{X), 

where U{X) := J7m(A) ••• Uo{X). Naturally the product ^/{^(xm; A) [/(A) depends 
analytically on A. 

Returning to the separate interval calculations on [£ . , .i:*], we see that by the 
procedure just outlined we can generate the solution approximations 

^±(a;*;A) = y±(±)(ar.;A) U^{X). 

The number of integration steps M can of course be different in each interval. 

Hence, after dropping the exponential prefactor and fixing the matching point to 
be a;*, the Evans function D{X;x^,) can be approximated by D{X;x*) where 

D{X;x^) := det (YL (x* ; A) Y+{x^-X)) 

= det(t/jo^(_j(a;*;A) yjt (+)(a;,; A)) • det [/"(A) • det C/+(A). 
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This is an analytic function of A. Indeed as we hinted previously, at each compu- 
tation step X = Xm we need only store j/jo^ [xm ; A) and the value 



n detC/,(A) 

l=m 

which gets updated at each step by simply multiplying the previous step value 
by the complex scalar determinental factor for current step. Hence to preserve 
analyticity for the Evans function using GGEM we must generate an approximate 
flow on Gr(n, k) ® C. 

6.5. Scaled GGEM. The scalar determinental factor just described, that we up- 
date at each step, grows exponentially. This would be tempered by the scalar 
exponential prefactor in the definition of the Evans function. An accurate practical 
procedure here is as follows (to be applied with due care). When integrating in the 
interval [^_,a::*], at each step, divide the scalar determinental factor in GGEM by 
exp((yLt|f (A)-l-- • ■ + ii'i^ {Xj)h) , where h is the stepsize, and the IJ,^ {)^) are the (spatial) 
eigenvalues, with positive real part, of A{—oo; A). When integrating in the interval 
[£+, x^,], at each step, divide the scalar factor by exp(— (^j*'(A) -I- • • • -I- A^^_fc(A))/i) , 
where the fJ,f{X) are the eigenvalues, with negative real part, of A(-\-oo; A). 

To see that this normalization is appropriate, we recall the Pliicker coordinates 
of Section [51 After applying the optimal Gaussian elimination algorithm to K^+i 
the i^^^th row elements of yi^^^ are themselves {n — k)k Pliicker coordinates; 
normalized by det Um+i- The ij^_|_]^th row elements of and det Um+l^ can be 

used to reconstruct the remaining Pliicker coordinates through the homogeneous, 
quadratic Pliicker relations. Hence the Pliicker coordinates, or complete set oi kx k 
minors, of Ym+i and differ by a factor det Um+i- It is well known that if the 

original vector field on V(n, k) is linear, then the Pliicker coordinates corresponding 
to Yjn+i satisfy a (larger) linear system of equations. In the left far field the 
coordinates thus grow exponentially, in fact with growth rate + ■ ■ ■ + /ijr(A); 

hence our recommendation to divide by the exponential factor suggested (with an 
analogous argument for the right far field). See Alexander, Gardner and Jones [2], 
Alexander and Sachs ^ , Brin [HI [17] or Allen and Bridges [5] for more details. 

7. Applications 

We present some numerical results for three different applications. The three 
applications reduce to the solution of a system showing multiple distinct exponential 
growth and decay rates in the stable and unstable subspaces, respectively. We show 
that our approach resolves this numerical obstacle successfully and can compete 
with the continuous orthogonalization method of Humpherys and Zumbrun j49j . 

7.1. Algorithms. We implement six different algorithms as follows. 

(1) Riccati-RK: Riccati method with fixed coordinatization with the flow of 
the Riccati vector field approximated by the classical fourth order Runge-Kutta 
method. We generically chose the coordinatization labelled by = {1, . . . , fc} and 
i+ = {fc -|- 1, . . . , n} for the left-hand and right-hand intervals, respectively. Hence 
if Y^{X) and Y^{X) denote the unstable and stable subspaces of A{—oo;X) and 
A{+oo; A), respectively, then we set 

y^{i±;X) = (ro±(A))(^±)„ .±(ro±(A))-'.±, 
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where denotes the i x i' submatrix of Y . We integrate the Riccati equation 

outhned in Corollary[l]in the two intervals and evaluate the modified Evans function 

^(A;..)^detf._/^- ^\(--^)y 

Provided neither Riccati flow becomes singular, this Evans function is analytic in 
the spectral parameter A. 

(2) Mobius-Magnus: Uses the Schiff and Shnider approach to integrate through 
singularities, combined with a Lie group method to advance the solution on the 
general linear group, as described at the end of Section [3] The same generic fixed 
coordinate charts are used as for the Riccati-RK method above. Over the inte- 
gration interval [xm, x^+i] with an equidistant mesh stepsize h, we advance the 
solution on the Lie algebra using the fourth order Magnus method 

= '^h{AixW)+Aixi^)) - 4h^[Aix^^lAix^^)], 

with the two Gauss-Legendre points (see Iserles, Marthinsen and N0rsett [55]) 

a;W = x„, + (5 - lVi)h and x\^} ^ x„, + (5 + ^V3)/i. 

We then compute the Mobius map ym+i = fJ-Hm ° ^^P to advance the solution in 
the fixed Grassmannian chart — for the left-hand interval = {I, . . . ,k} while for 
the right-hand interval i+ = {fc -I- 1, . . . , n}. We evaluate the same Evans function 
as for the Riccati-RK method above. 

(3) GGEM-RK: Scaled Grassmann Gaussian elimination method, with the clas- 
sical fourth order Runge-Kutta method used to advance the solution on the Stiefel 
manifold, as described in Sections [5] and 16.51 We evaluate the Evans function 
D{\\x*) in Sectionini 

(4) GGEM-LG: Same as GGEM-RK but with a fourth order Magnus method 
used to advance the solution on the Stiefel manifold instead, i.e. Ym+i = exp((Tm) 
where am is generated as for the Mobius-Magnus method above. 

(5) Riccati- QOGE: Riccati method with coordinate swapping as described in 
Section 321 We have chosen to implement the method in the following form. At 
each integration step we advance the solution on the Stiefel manifold using the Mag- 
nus method (we could also use a Runge-Kutta method here). We apply elementary 
column operations to the resulting solution matrix to convert the pre-determined 
rows indexed by i from the previous step to the identity matrix. Then if ||y||oo is 
less than or equal to a tolerance size, we keep this index t for the next step. If it 
is greater, we apply QOGE at the end of the next step after advancing the solu- 
tion on the Stiefel manifold, thus generating a new index. As for GGEM-RK and 
GGEM-LG, we update the scalar determinental factor at each step (produced by 
the elementary column operations with the pre-determined index or QOGE). We 
divide the scalar determinental factor by the scalar exponential factors, as described 
for the scaled GGEM method. We evaluate the same Evans function also. 

(6) CO-RK: Continuous orthogonalization method of Humpherys and Zumbrun 
with the classical fourth order Runge-Kutta method used to advance the solu- 
tion on the Stiefel manifold of orthonormal frames. The initial conditions Qq{X) 
for the Q- matrices are obtained by QR-factorization of Y^{X). From Humpherys 
and Zumbrun [49j . to ensure analyticity we must also solve the scalar problems 
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Figure 1. The Evans function of the Boussinesq system for the 
unstable pulse having wave speed c = 0.4. The left plot is gen- 
erated by the Riccati-RK method with fixed coordinate patches 
identified by = {1,2} over [-8,0] and i+ = {3,4} over [0,8]. 
The right plot shows the result of the CO-RK method. 

(deti?±)' = Tr(QtA(a;;A)Q-(g±(A))U(±(X); A)Qj(A)a;) deti?±. The Evans func- 
tion is then given by 

i:i(A;a;*) = deti?-(x*; A) • det A) • det(Q-(a;*; A) Q+{x^;X)). 

7.2. Boussinesq system. As the first test system, we consider the Boussinesq 
system studied by Humpherys and Zumbrun ^49]. The (good) Boussinesq equation, 
expressed in a co-moving frame moving to the right with wave speed c, is given by 

Utt = (1 - C^) Uxx + '2cUxt - Uxxxx - {U^)xx- 

It has solitary wave solutions of the form u{x) = — c^)sech^(i-\/l — x), where 
|c| < 1. These waves are stable when 1/2 < |c| < 1 and unstable when |c| < 1/2. 

If we consider small perturbations about the travelling wave u we generate a 
linear spectral problem of the form Y' — A{x\ A)y, where 

^0 1 0^ 

10 

1 

\-A2-2m" 2Ac-4m' (1-c2)-2u 0/ 

When the spectral parameter A lies in the right-half complex plane the eigenvalues 
of A(±cx); A) spectrally separate into two growth and two decay modes, i.e. k = 2. 
We used £± — ±8 in our experiments. 

In Figure [T] we show the Evans function computed along the real axis from 
A = to A = 0.2 for the unstable pulse with c = 0.4. The Riccati-RK (left plot) 
and CO-RK (right plot) methods detect a zero of the Evans function near A = 
0.155, indicating an unstable eigenvalue there. An accurate value of the eigenvalue 
can be found by using a standard root-finding method. This yields the values in 
Table [TJ The Riccati-RK and the CO-RK methods both converge to the same 
eigenvalue when the number of steps N increases. As a check, the Matlab ode45 
solver was used with a relative tolerance 10~* and absolute tolerance 10"^" to 
integrate the systems, leading to the same resulting eigenvalue: A = 0.15543141 for 
both methods. 



A{x;X) = 
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Table 1. Zero of the Evans function for the Boussinesq prob- 
lem computed with the Riccati-RK method (with fixed coordinate 
patches identified by i =~ {1,2} over [—8,0] and i"*" = {3,4} over 
[0,8]) and the CO-RK method. Here iV is the number of (equidis- 
tant) steps used in the mesh. 



N 


Riccati-RK 


CO-RK 


128 


0.15544090 


0.15540090 


256 


0.15543184 


0.15542952 


512 


0.15543143 


0.15543129 


1024 


0.15543141 


0.15543140 


2048 


0.15543141 


0.15543141 


4096 


0.15543141 


0.15543141 
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Figure 2. Error in the eigenvalue, vs stepsize (upper panel) and 
vs cputime (lower panel), matching at = 0. 



Function evaluation for the Riccati vector field requires three matrix-matrix 
multiplications. This is the same number of matrix-matrix multiplications needed 
to evaluate the Drury-Oja vector field. However, the matrices in the Drury-Oja 
vector field are nx k and nx {n — k), respectively, while the matrices in the Riccati 
vector fields have smaller dimension: {n — k)x k. Because of the smaller dimension 
of the systems to be integrated, our Riccati approach is faster than the continuous 
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Figure 3. Error in the eigenvalue vs stepsize (upper panel) and 
vs cputime (lower panel), matching at = 8. 



orthogonalization problem. For example, to construct Figure [T] the Evans function 
was evaluated at 200 distinct A values between A = and A = 0.2. Using the 
fourth-order Runge-Kutta method with TV = 512 steps, this required 33 seconds 
for the CO-RK method, while the Riccati-RK method needed 24 seconds (Matlab- 
implementation, CPU 2.4GHz). 

In Figure [2] we compare the error in the eigenvalue and efficiency of computation 
for all six methods, when we match at x, = 0. We see that the methods that use 
the Magnus expansion to advance the solution on the Stiefel manifold are the most 
accurate for a given stepsize. They are also the most efficient, delivering the best 
accuracy for given computational effort. The Riccati-RK method does not suffer 
from singularities for the chosen fixed patches when matching at x* = 0, at least 
for the range of values of the spectral parameter in the vicinity of the eigenvalue (as 
well as the origin and anywhere inbetween). However, if we change the matching 
point to X* = = +8 there are singularities in the Riccati-RK solution (as a result 
of a vanishing determinant ofu~). In particular, a singularity appears around x = 2 
for A equal to the eigenvalue (see Figure [3]) . Hence we compare the remaining five 
methods in Figure [3] in this case. We see that using the Mobius-Magnus method to 
integrate through a singularity does not introduce loss in accuracy. The resulting 
Evans function can have poles, as seen in Figure [H which appear at A- values where 
the Riccati equation has a singularity at the matching point. This means that in 
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Figure 4. The Evans function when the Mobius-Magnus method 
is apphed (left panel) with the matching point as a;* = +8. The 
entries of the solution , passing through the singularity when 
integrating from —8 to 8, for A = 0.15543141 (right panel). 




Figure 5. The Evans function which results when the GGEM- 
RK scheme is applied over [—8,8], matching point in a;* = 8 (left 
panel). The entries of for A = 0.15543141 (right panel). 




Figure 6. The Evans function when the Riccati-QOGE method 
is applied (left panel). The entries of yio are shown for A = 
0.15543141 (right panel). The criterion used for swapping to a 
new coordinate patch was ||j/^||oo > 2. 
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Figure 7. Error in the eigenvalue for different choices of the 
matching point. The number of steps in the equidistant mesh was 
N = 512. 
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Figure 8. |-D(A)| for A equal to the eigenvalue for different match- 
ing points. The number of steps in the equidistant mesh was N — 
512. 



some cases the matching point should be chosen rather carefully in order not to 
have the poles interfering with the eigenvalue(s). When applying GGEM-RK the 
Evans function is analytic and the choice of the matching point is less important. 

Figure [5] shows the Evans function obtained when the GGEM-RK evolves yi from 
£- — —8 to = -|-8. To construct the plot in Figurc[5l the quasi-optimal Gaussian 
elimination process was applied at each step in the integration. However it is clear 
from the right plot in Figure El that multiple successive steps can be integrated 
in the same coordinate patch. For example, between x — —8 and x — —3 the 
coordinate patch does not change. Performing the whole quasi-optimal Gaussian 
elimination process only when a certain criterion is satisfied, reduces the comput- 
ing time. Using the Riccati-QOGE method, we change the coordinatization when 
||y~||oo > 2. This generates an Evans function very similar to that in Figure [5l As 
seen in Figure [6] the quasi-optimal Gaussian elimination process is then performed 
only two times for A equal to the eigenvalue. 
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Figure 9. The travelling wave solution for (5 = 0.1 and m = 9. 



We compare the error in computing the eigenvalue for different choices of match- 
ing point — in fact for anywhere in the interval [—8, 8] — for all six methods in 
Figure [7l We see that the most accurate and robust methods are the GGEM-LG 
and Riccati-QOGE methods. Some methods, such as the Riccati-RK method as 
discussed already, break down when singularities impinge on the matching point — 
the singularities in the Evans function are observed in Figure ID Generally we also 
see in Figure [7] that the GGEM-LG and Riccati-QOGE methods outperform the 
CO-RK method in terms of accuracy. 

Overall, we observe in this example that when computing the eigenvalue, those 
methods based on the Magnus expansion are superior in accuracy and efhciency. 
Note that for GGEM-RK, the quasi-optimal Gaussian elimination process is an 
additional nfc^ operation. However, to ensure analyticity for the CO-RK method, 
there are two additional matrix-matrix multiplications in the equations for det 
(operational cost kn^) required at each step. 

7.3. Autocatalytic fronts. As a second example, we study travelling waves in a 
model of autocatalysis in an infinitely extended medium 



Here u{x,t) is the concentration of the reactant and v{x,t) is the concentration 
of the autocatalyst. We suppose {u, v) approaches the stable homogeneous steady 
state (0,1) as X ^ — oo, and the unstable homogeneous steady state (1,0) as 
X — + +00. The diffusion parameter 6 is the ratio of the diffusivity of the reactant 
to that of the autocatalyst and m is the order of the autocatalytic reaction. The 
speed of the co-moving reference frame is c. The system is globally well-posed for 
smooth initial data and any finite 5 > and m > 1. 

From Billingham and Needham |9j we know that a unique heteroclinic connection 
between the unstable and stable homogeneous steady states exists for wavespeeds 
c > Cmin. The unique travelling wave for c = Cmin converges exponentially to the 
homogeneous steady states and is computed by a simple shooting algorithm (see 
Balmforth, Craster and Malham [8|). The resulting travelling wave for 5 = 0.1 and 
m = 9 is shown in Figure [H 



ut = 5uxx + cux — uv' 



rn 



Vt = Vxx + CVx + UV 



m 
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Figure 10. Zero contour lines of the real (solid) and imaginary 
(dashed) parts of the Evans function for the autocatalysis problem 
with 5 = 0.1 and to = 8 (left panel) and m = 9 (right panel). 
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Figure 11. Error in the eigenvalue in the first quadrant when 
6 — 0.1 and m = 9 for different methods, vs stepsize (upper panel) 
and vs cputime (lower panel), matching at = 0. 



The stability of the travelling wave of velocity c can be deduced from the location 
of the spectrum of the eigenvalue problem Y' = A{x; X)Y , where 

1 ^ \ 

1 

muv"'~^/S -c/6 ' 

A — TOuw™^^ —cj 

where u and v represent the travelling wave solution. 

The pulsating instability occurs when i5 < 1 is sufficiently small and to is suf- 
ficiently large (see Metcalf, Merkin and Scott [72] and Balmforth, Craster and 



A{x;X) = 






X/5 + v"'/6 
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Figure 12. Error in the eigenvalue in the first quadrant when 
i5 = 0.1 and m — 9, for different methods, for different matching 
points. The number of steps in the equidistant mesh is iV = 256. 



Malham [5]). For 6 fixed and m increasing, a complex conjugate pair of eigen- 
values crosses into the right-half A-plane signifying the onset of instability via a 
Hopf bifurcation. Figure [TO] shows the onset of this instability for 5 = 0.1 as m is 
increased from 8 to 9 (see Aparicio, Malham and Oliver [7]). The figure shows the 
zero contour lines of the real and imaginary parts of the Evans function. Solid lines 
correspond to zero contours of the real part of D{X), dashed lines to the imaginary 
part of D{X). We see that a complex-conjugate pair of eigenvalues crosses into the 
right-half plane, indicating the onset of instability. Figure fTOl was constructed using 
the Riccati-RK method with the fixed coordinate patches identified by ={1,2} 
from X = —10 to x^, = —7, and i"*" — {3, 4} from x = 10 to = —7. The matching 
point X* = — 7 is chosen roughly centred on the wavefront. 

We compare the Riccati-RK, Mobius-Magnus, CO-RK and GGEM-LG methods 
in Figure[TT]where we plot the absolute error in the eigenvalue vs the stepsize (upper 
panel) and also vs cputime (lower panel) . The eigenvalue in question is that in the 
first quadrant in Figure [TO] for 5 = 0.1 and m = 9. Figure [11] was generated 
as follows. Starting with an initial guess lying within a small square around the 
eigenvalue, we iterated a standard root finding algorithm until we arrived in a 
square (containing the eigenvalue) which was smaller than a preset tolerance. We 
see in Figure [TT] that the Riccati-RK method produces a slightly better error for a 
given stepsize, and is marginally more efficient than the GGEM-LG method. The 
CO-RK method produces a larger error for a given computational effort. This is 
not surprising, as again, the matrices in the Drury-Oja vector field are twice as big 
(4 X 2) as the ones in the Riccati vector fields (2x2). 

When we match at = — 7, there is little to distinguish the Riccati-RK, Mobius- 
Magnus, CO-RK and GGEM-LG methods. We compare all four methods for dif- 
ferent matching positions x* in Figure [T^J which was generated using the same 
root finding criteria as for Figure [TTl except all the methods used N = 256 steps. 
Note that the errors in the Mobius-Magnus and GGEM-LG methods are uniform 
for any matching values in the range [—10, 10]. The CO-RK method error doesn't 
vary that much either and is slightly larger. Note that no values are plotted for 
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Figure 13. Contour lines of ji'lA)! for 8 — 0.1 and m = 9 when 
using the Riccati-RK, Mobius-Magnus, CO-RK and GGEM-LG 
methods (order: top three down to bottom three), matching at 
positions cc* = —8,0, +8 (left to right). 



the CO-RK method at the matching points = 8, 10. In these cases the clas- 
sical Runge-Kutta method applied to the Drury-Oja vector field is unstable for 
N — 256 steps for some A-values close to the eigenvalue. This problem is resolved 
by increasing the number of steps to N — 512. For a range of matching positions 
roughly in [—10, 0], there are no singularities of the Riccati-RK method in the left 
and right-hand integration intervals for values of the spectral parameter A close to 
the eigenvalue. Indeed for this range of matching positions the Riccati-RK method 
delivers the best accuracy. However for matching positions outside this range, for 
values of the spectral parameter A not far from the eigenvalue, the Riccati-RK so- 
lution does have a singularity for some matching points (which we can see in the 
contour plots in Figure [T3|) . This makes the eigenvalue-searching algorithm fail — 
indicated by no error points for those matching position values. We also do not 
show the error for the Riccati-RK method for the matching points = —8, —10, as 
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□ 



Figure 14. We show three different closed contours in the first 
quadrant of the complex A-plane in each of the left panels. In the 
top two left panels, the contour encloses the eigenvalue in that 
quadrant (found in Figure [T51 for S = 0.1, m = 9). In the corre- 
sponding three panels on the right we show how the argument (in 
multiples of 2tt) of the Evans function D{X) changes as we perform 
a complete circuit of the contour. The Evans function was com- 
puted using the GGEM-LG method matching at x^, = +14. The 
number of patch changes performed for each fixed A value are in- 
dicated. 



there are singularities in the Evans function close to the real axis for these match- 
ing points (again see Figure [T3)) . This means that we cannot for example, apply 
the argument principle in the first quadrant, though starting sufficiently close to 
the eigenvalue we can still use the Riccati-RK method as part of a root-finding 
algorithm to determine the eigenvalue. 

Figure [Ql shows the contour lines of |-D(A)| for 6 = 0.1 and m = 9, close to 
the eigenvalue in the first quadrant, when using the Riccati-RK, Mobius-Magnus, 
CO-RK and GGEM-LG methods, respectively, and matching at three positions 

= -8, 0, +8. We see that the CO-RK and GGEM-LG methods show the least 
sensitivity to the choice of matching position and produce smooth contour plots 
for all three matching points. The contour plots across shape and scale look very 
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similar for both these methods. By contrast the Riccati-RK and Mobius-Magnus 
methods appear to develop singularities close to the eigenvalue when the matching 
position X* is or +8. 

Lastly in Figure [14] we demonstrate the argument principle for counting zeros 
of the Evans function inside closed contours. We computed the Evans function 
using the GGEM-LG method and matched at a;* = +14. As expected, if the closed 
contour in the complex A-plane encloses the eigenvalue, then the change in the 
argument of the Evans function around the complete contour is one (once we have 
accounted for the 2tt factor in the argument principle). We also show, for each fixed 
A value, the number of patch changes that occured as we integrated from a; = — 14 
through to a; = +14. 



7.4. Ekman boundary layer. The third test system is a boundary layer flow over 
a flat plate which is infinitely extended in the x and y direction and rotates around 
the half infinite z-axis with a given rotational speed. Linear stability of the Ekman 
boundary layer has been investigated in Allen and Bridges [}B| and Allen [4J using 
the compound matrix method. The fiow is governed by the continuity equation 
Ux + Vy + w z =0, and the Navier-Stokes equations in a co-rotating frame 



Ut + UUx + VUy + WUz + -^Px 
Vt + UVx + VVy + WVz + -^Py ■ 
Wt + UWx + VWy + WWz + 



Ro 
2 



RoE, 



-Pz 



R„ ("2:2; - 
^{Vxx - 
^{Wxx 



' ^ RT' 



Here Ro, Ro and Ek denote the Reynolds, Rossby and Ekman numbers, respectively. 

After non-dimensionalization and setting Rg = Ro, E^, = 1, the linear stability 
of the boundary layer is determined by the eigenvalues A of the linear problem 
Y' ^ A{z; X)Y, where (see Allen [3 p. 176]) 



Aiz;X) 



( 



— a(z. A) 






1 



6(z, A) 











6(z,A)-72 






-2 
1 



and 



V{z) = cos(e)(l — exp(— z) cos(z)) + sin(e) exp(— z) sin(z), 
Vz (z) — exp(— z) (sin(z + e) + cos(z + e)) , 

C/(z) = — sin(e) (l — exp(— z) cos(z)) + cos(e) exp(— z) sin(z) , 

Vzz{z) — —2 exp(— z) cos(z + e), 

a(z. A) = 7^ + iRc7'(7C^(z) - + i7Rc[/..(z), 

6(z,A) = 272 + R,(z7C/(z) + A). 

Here the parameters 7 and e represent the radial and angle components, respec- 
tively, of a polar coordinate parameterization of horizontal wavenumbers associated 
with the X and y directions — see Allen and Bridges [6] for more details. 

We choose the fixed coordinate patch identified by i+ = {1, 2, 3} to integrate the 
corresponding Riccati equation from z: = £+ to z* =0. The boundary condition for 
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Figure 15. Neutral curves for rigid wall, Ro fixed (values indicated). 





the rigid wall at z* = as given in Allen and Bridges is 

yi(0;A) = y2(0;A) =r5(0;A) = 0. 

If 

/I 

= 1 

\0 
then the boundary conditions are equivalent to 

dct(r- • r+(0; A)) = 0. 

We thus compute the Evans function: 

DiX;z,) = det(Y- ■ f{\\ = dot ( 1 I = y23(^*; A). 
V ^^^/ \y2i y22 hi I 

We computed neutral curves, i.e. curves in the e-7 plane where i?e(A) = 0, 
using the Riccati-RK method with i+ = {1,2,3} to compute y(0; A) € C^^'^ and 
consequently the Evans function -D(A; z^^,). For continuation of the curves we used 
the Matlab package MatCont which uses pseudo-arclength continuation (Dhooge, 
Govaerts and Kuznetsov [29]). Figures fT5l and [T6l show the neutral curves which 
match those in Allen and Bridges and Allen [4] . The integration of the Riccati 
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Figure 17. Contour plots of |D(A; 0)| for Re = 140, e = 0.014156, 
7 — 0.70575. An equidistant mesh was used with 500 intervals and 
values computed on a 20 x 20 grid in the complex A-plane. 



system has been done with the Matlab ODE-solver ode23s from z = 10 to = 
(as in Allen and Bridges) with absolute and relative tolerances 10~^ and 10^^. The 
stable subspace of A{+oo;X) was constructed using the Matlab eigenvalue-solver 
eig (we also used direct formulae for the eigenvectors to construct analytic bases for 
the stable subspace but this did not significantly change the overall performance). 

For comparison we also implemented the CO-RK and GGEM-LG methods. We 
tested the performance of all three methods, in each case evaluating the Evans 
function on a 20 x 20 grid for A in the complex plane. In Figures [17] we present 
contour plots of \D{X; 0)|, and see that all methods find a root at A « 0.002 — 0.1171. 
The computation times for a 2.4Ghz machine were: 134 seconds for Riccati-RK, 
150 seconds for GGEM-LG and 189 seconds for CO-RK. As a comprehensive check, 
we also implemented the compound matrix method (i.e. Pliicker coordinates, of 
which there are 20), described in Allen and Bridges, for this performance test. As 
expected, since this method involves integrating a linear system of order 20, it was 
an order of magnitude slower (while giving the same results). 

8. Concluding discussion 

We have shown that the new scaled Grassmann Gaussian elimination method as 
well as the Riccati method with quasi-optimal patch swapping, compete with the 
continuous orthogonalization method for computing the Evans function. Both new 
methods deliver superior accuracy for the same computational cost when combined 
with Lie group Magnus integration to advance the solution. Moreover, as hoped, 
numerically these new methods appear to be robust in the sense that they are 
insensitive to the choice of the matching position in the computational domain. We 
now outline several directions in which we plan to use and extend these methods. 

One of the main goals we have had in mind in this paper is that of large scale spec- 
tral problems, in particular the stability of travelling waves with a multi-dimensional 
structure. There is recent research extending the Evans function approach in 
this direction — see Deng and Nil ^7\, Gesztesy, Latushkin and Makarov [3^ and 
Gesztesy, Latushkin and Zumbrun |37| . From a numerical perspective we have, 
together with Niesen, implemented some of the methods we propose in this paper 
in a multi-dimensional context. In particular, it is well known in autocatalysis and 
combustion that planar travelling fronts can be unstable to transverse perturbations 
and develop into steadily propagating travelling fronts with wrinkles. In Ledoux, 
Malham, Niesen and Thiimmler [66j we show that the wrinkled fronts themselves 
develop an instability as a diffusion parameter is further increased. 
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For large scale problems the Lie group methods we propose using the Mag- 
nus expansion may become prohibitive. This is because of the effort required to 
compute the matrix exponential — see Moler and Van Loan , Celledoni and Iser- 
les , Munthe-Kaas and Zanna 78] and Iserles and Zanna [53, • For the examples 
we considered this was not an issue. However it remains to be seen if such Lie 
group methods will be cost effective for larger problems — the methods we proposed 
based on Runge-Kutta integration such as GGEM-RK can be used as they scale 
favourably with system size. 

The constructs and Grassmannian reductions we have considered in this paper, it 
turns out, have their origins in the control theory literature dating back to the early 
seventies, in particular in the pioneering papers of Hermann and Martin O |45l 
ESI mi [48], Martin and Hermann [71] and Brockett and Byrnes [18]. We also found 
Bittanti, Laub and Willems [12], Lafortune and Winternitz [G^, Rosenthal 87], 
Shayman [jQl and Zelikin |100| particularly useful resources. A future direction 
we would like to explore is whether there are any applications of the numerical 
methods we have outlined here to practical non-autonomous control problems? 

Riccati methods in particular also have their origins in the quantum chemistry 
literature also dating back to the early seventies — a recent survey of these numer- 
ical methods can be found in Chou and Wyatt [21]. However also see Light and 
Walker [S7], Johnson [S5], Hutson [ST] and Gray and Manopoulous 139]. In par- 
ticular the log-derivative and i?-propagation methods correspond to special choices 
of Grassmannian patch in the Riccati methods we mention above. Priifcr meth- 
ods, for which we can think of the patch evolving, originate even further back; see 
Priifer [84, and Pryce ^85]. 

Of course, our quasi-optimal Gaussian elimination process for choosing a suit- 
able representative patch was inspired by the Schubert cell decomposition of the 
Grassmann manifold; see for example Billey flW , Griffiths and Harris [40] , Kleiman 
and Laksov 58j, Kresch [60], Postnikov [83 , Sottile [93^ and, in a somewhat dif- 
ferent vein, Kodama [59] . Since the Grassmann manifold is the disjoint union of 
Schubert cells, the question is, can we express the flow on the Grassmann manifold 
as a flow on Schubert cells (see Griffiths and Harris and also Ravi, Rosenthal and 
Wang [86])? Can we construct the corresponding flow on the cohomological ring of 
Schubert cycles (Chern 22]; Fuhon [35])? 
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